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Abstract 

The Quantum Unique Ergodicity (QUE) conjecture of Rudnick- 
Sarnak is that every eigenfunction /> n of the Laplacian on a manifold 
with uniformly-hyperbolic geodesic flow becomes equidistributed in 
the semiclassical limit (eigenvalue E n —> oo), that is, ‘strong scars’ 
are absent. We study numerically the rate of equidistribution for a 
uniformly-hyperbolic Sinai-type planar Euclidean billiard with Dirich- 
let boundary condition (the ‘drum problem’) at unprecedented high 
E and statistical accuracy, via the matrix elements {(j) n , Acf> m } of a 
piecewise-constant test function A. By collecting 30000 diagonal ele¬ 
ments (up to level n ~ 7 x 10 5 ) we find that their variance decays with 
eigenvalue as a power 0.48=1=0.01, close to the estimate 1/2 of Feingold- 
Peres (FP). This contrasts the results of existing studies, which have 
been limited to E n a factor 10 2 smaller. We find strong evidence for 
QUE in this system. We also compare off-diagonal variance, as a func¬ 
tion of distance from the diagonal, against FP at the highest accuracy 
(0.7%) thus far in any chaotic system. We outline the efficient scaling 
method used to calculate eigenfunctions. 

*now at Department of Mathematics, Dartmouth College, 6188 Bradley Hall, Hanover, 
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1 Introduction 


The nature of the quantum (wave) mechanics of Hamiltonian systems whose 
classical counterparts are chaotic has been of long-standing interest, dating 
back to Einstein in 1917 (see for a historical account). The field now 
called ‘quantum chaos’ is the study of such quantized systems in the short 
wavelength (semiclassical, h —► 0 or high energy) limit, and has become 
a fruitful area of enquiry for both physicists (for reviews see |22, [3D]) and 
mathematicians (see |4311441 551 1 in recent decades. In contrast to those of 
integrable classical systems, eigenfunctions are irregular. A central issue, and 
the topic of this numerical study, is their behavior in the semiclassical limit. 

We consider billiards, a paradigm problem in this field. A point particle 
is trapped inside a bounded planar domain flcK 2 and bounces elastically 
off the boundary T = dfl. Its phase space coordinate is (r, 6) e x S 1 , 
with position r := (x,y) and (momentum) direction 9. The corresponding 
quantum-mechanical system is the spectral problem of the Laplacian in 
with homogeneous local boundary conditions which we may (and will) take 
to be Dirichlet, 

E n (j) n , (1) 

0n(r) = 0. (2) 

Eigenfunctions 0 n are real-valued and normalized 

(07i) 0m) • / 0n(r)0m(r)dr — $nr?i) (3) 

Jn 

where dr := dxdy is the usual area element, and the corresponding ‘energy’ 
(or frequency) eigenvalues are ordered E\ < E 2 < E 3 < ■ ■ ■ oo. We will also 
write Ej = k | where kj is the wavenumber. This Dirichlet eigenproblem [ 32 ], 
also known as the membrane or drum problem, has a rich 150-year history of 
applications to acoustics, electromagnetism, optics, vibrations, and quantum 
mechanics. 

When the classical dynamics (flow) is ergodic it is well-known that almost 
all eigenfunctions are ‘quantum ergodic’, in the following sense. 

Theorem 1 (Quantum Ergodicity Theorem (QET) |46L 1551121L162] 1 

Let fl e M 2 be a 2D compact domain with piecewise smooth boundary whose 
classical flow is ergodic. Then for all n except a subsequence of vanishing 
density, 

(0 n , A<f n ) - A —> 0 asn ^ oo, (4) 
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for all well-behaved functions A : > M. 

The operator A is multiplication by A, a ‘test function’ whose spatial average 
is A := vol ^ f n H(r) dr. The physical interpretation is that almost all quan¬ 
tum expectation values (diagonal matrix elements) of the observable A tend 
to their classical expectation A, an example of the Correspondence Princi¬ 
ple of quantum mechanics. Note the choice ‘almost all’ need not depend on 
A gSJ. Equivalently, almost all probability densities \4>j\ 2 tend to the uni¬ 
form function , weakly in L 1 (f2) (see Fig. 0). The proof relies on the 
machinery of semiclassical analysis including Fourier Integral Operators (see 
[EBUM! and references within). In our context the ‘well-behaved’ requirement 
is, loosely speaking, that A not be oscillatory on the (vanishing) wavelength 
scale 1/kj. (Formally A must be a zeroth-order pseudo-differential operator, 
a Weyl quantization of the principal symbol H(r)). We investigate only mul¬ 
tiplication operators, that is, no dependence on momentum, and, as we will 
see below, by further restricting to piecewise constant test functions we will 
exploit a huge numerical efficiency gain. 

Our study is motivated by the fact that QET tells us nothing about the 
rate of convergence of (cf) n , A<f> n ) or the density of the excluded subsequence, 
both of which are needed to understand the practical applicability of the 
QET in quantum or other eigenmode systems. We are interested in how the 
size of the deviation (0 n , A(p n ) — A varies with eigenvalue E n . We define its 
‘local variance’ (mean square value) at energy E by 


V a (E) : = 


1 

ME) 


£ 

n:E n e[E,E+L(E)} 


2 


{4>n, A<f n ) - A 


( 5 ) 


where Nl(E) N(E + L(E)) — N(E) and the level counting function is 


N(E) := #{n : E n < E} (6) 

Here we envisage an energy window width L(E) = 0(E F 2 ), that is, a 
wavenumber window of width 0(1): this contains 0(E 1 / 2 ) eigenvalues by 
Weyl’s Law For practical reasons (Section 14.11) we will in fact use 

other L(E) widths which nevertheless contain many (« 10 3 ) eigenvalues. We 
will test the following asymptotic form for the variance. 

Conjecture 1 (Power-law diagonal variance) For ergodic flow, as E —> 
oo, there is the asymptotic form, for some a and 7 , 


V A (E) ~ aE~\ 


( 7 ) 
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A random-wave model for eigenfunctions (Section 12. Ill predicts 7 = 1/2 and 
a certain prefactor a RW . We will also test a more elaborate heuristic from 
the physics literature (Section 12.211 which predicts 7 = 1/2 and the following 
different prefactor. 

Conjecture 2 (Feingold-Peres diagonal variance [26,123]) For ergodic 
flow with no symmetries other than time-reversal, 

v a(E) ~ ^^ W 172 ( 8 ) 

vol(S2) 

where the symmetry factor is g = 2. 

Here the prefactor Ca(uj) is the power spectral density (also known as clas¬ 
sical variance |44j . spectral measure EH ES], or fluctuations intensity |20|). 


/ OO 

C A {r)e iu}T dr, (9) 

-OO 

the time autocorrelation of A being 

Ca(t) := lim ^ f A(r(t))A(r(t + r)) dt, (10) 

T^oo 1 J Q 

where r(i) is any uniformly-distributed (ergodic) trajectory. 


Remark 1.1 In the physics literature much stronger conjectures are often 
discussed, such as individual matrix elements (cf) n , A<p n ) being pseudo-randomly 
distributed with variance given as above. We present evidence for this in 
Section \j\ At the other extreme, proven theorems (and some numerical 

work w often involve sums of the form S P (E ; A) := N(E) 1 JA-e,<e Ucp^Acpfl 

Note that asymptotic decay S P (E ; A) ~ bE ~ p7 / 2 is equivalent to Conjecture Q 
However we will study Va(E) rather than S P (E ; A) for the following impor¬ 
tant practical reasons: 


A 


1. Narrow spectral windows are needed rather than the complete spectrum, 
allowing much higher eigenvalues to be included in the statistics. 

2. Asymptotic behavior will emerge sooner since data at high E n are not 
averaged with that from the lower part of the spectrum. 
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Figure 1: Billiard and test function used in this study, a) Quarter generalized 
Sinai billiard formed by two circular arcs which meet at (1,1) at angles 
61 = 0.4 to the horizontal and 6 b = 0.7 to the vertical. The straight sections 
lie on the axes and meet the arcs at right angles, b) Piecewise-constant A(r) 
which takes the value 1 inside the region and zero elsewhere. 8 Qa \ T is 
a straight line which does not meet either wall at right angles. 


About 10 % agreement 1 with Conjecture |21 has been shown in the quantum 
bakers map EH and rough agreement with 7 = 1/2 has been found in 
hyperbolic polygons [2]. For general ergodic flow, proven bounds on the 
power-law rate 7 in Conjecture |T| are quite wide: Zelditch J3H] has shown 
that S P (E;A ) = 0((logl?) _P//2 ), implying 7 > 0, and shown that for 
generic A (more precisely, one with nonzero mean sub-principal symbol), 
7 < 1 . 

There are special ‘arithmetic’ manifolds with ergodic flow for which pow¬ 
erful number-theoretic tools pEil, fPJ allow much more to be proven. For 
the quotient manifold H 2 /AL(2,Z) Luo and Sarnak ^ recently showed 
S 2 {E-A) r-u B{A)E 1//2 , where the prefactor is a quadratic form B(-) di¬ 
agonalized by the eigenfunctions themselves. It takes the value B(^> n ) = 

1 These researchers also studied the hydrogen atom in a strong magnetic field (nearly 

completely ergodic with sticky islands), and found some agreement at the 20% level, 
however they admit that the agreement was ‘unexpectedly good’ since it depended on a 
choice of smoothing parameter. 
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Figure 2: Illustration of spatial equidistribution of \4> n \ 2 (shown as density plots; 
larger values are darker) with increasing eigenvalue: n th Dirichlet eigenfunction 
for mode numbers n = 1,10,100,1000 and n ~ 50000. 

C(0)L(i, (p n ), where L(-, •) is an L-function |il] . Thus the power law appear¬ 
ing in Conjecture El holds but the prefactor differs from g — 2. This is hardly 
surprising; Sarnak (Elj notes that a simple reflection symmetry is enough to 
cause g = 0. Arithmetic systems are very special and have many symme¬ 
tries: all periodic orbits possess a single Lyapunov exponent, and eigenvalue 
spacing statistics are unusual for ergodic systems DEI- This makes the study 
of a planar Euclidean billiard system without symmetry (where Lyapunov 
exponents differ), for which no number-theoretic analytic tools exist, partic¬ 
ularly interesting. We choose such a generic billiard (with hyperbolic, i.e. 
Anosov, flow [17] ) for our numerical experiments (see Fig. |TJ) . Numerical 
tests of Conjecture |T| (in the form of S\(E] A), see Remark ll.il) exist for low 
eigenvalues (n < 6000) in the Anosov cardioid billiard [5J: various powers 
were found in the range 7 = 0.37 to 0.5, and up to 20% deviations from the 
prefactor in Conjecture El 

We address two other questions, the first regarding the excluded sub¬ 
sequence in Theorem El the second regarding off-diagonal matrix elements. 


Conjecture 3 (Quantum Unique Ergodicity (QUE) [?TJ) There is no 
excluded subsequence in Thm. |7] 

‘Unique’ refers to the existence of only one ‘quantum limit’ (any measure 
to which | (f ) n |“ tends weakly). Made in the context of general negatively- 
curved manifolds, this conjecture has remarkably been proved for arithmetic 
manifolds m- In contrast, QUE has been proven not to hold for certain¬ 
dimensional quantizations of Arnold’s cat map (23 ■ This begs the question: 
for which ergodic billiards, if any, does this conjecture hold? For instance 
there is strong evidence mm (^t no proof [23] ) that a sequence of ‘bouncing 
ball’ modes (eigenfunctions concentrated on cylindrical or neutrally stable 
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orbits which occupy zero measure in phase space) persists to arbitrarily high 
eigenvalues in ergodic billiards such as Bunimovich’s stadium. Since such 
modes are not spatially uniform, the conjecture would not hold for this shape. 
There are also more subtle nonuniformity effects: in the physics community 
enhancements of \(j) n \ 2 , dubbed ‘scars’ by Heller 221 GUI], are known to exist 
along isolated unstable periodic orbits, and there has been a long-standing 
debate as to whether these may prevent QUE. In Section 14.41 we discuss how 
our results relate to scarring. 

Finally we consider the size of off-diagonal matrix elements (<j) n , A(ft m ). 
Define A k{E) := 2n/[E 1 t 2 vol(D)] as the mean level spacing (in wavenumber). 


Conjecture 4 (Feingold-Peres off-diagonal variance [2B, iQJ) Fixuj 6 
M. Then for ergodic flow, as E —> oo there is the asymptotic result, 


Va(E-u) 


A k (E) 

2 e(E)N L (E) 


^2 |(0n, Mm)? 

m,n:En£[E,E-\-L(E)] 

| km k n uj | ^e(£/) 


r\j 


where 0 < e(E) = 0(E~ 1 ^ 4 ). 


^Hn-l/2 

vol(fl) 


( 11 ) 


Va(E\uj) measures the mean off-diagonal quantum variance a ‘distance’ (in 
wavenumber units) u from the diagonal, which is thus given by classical vari¬ 
ance at frequency u (see Section 12.21 for a heuristic argument). As above, the 
choice L(E ) = 0(E 1//2 ) is envisaged. The rate at which the wavenumber win¬ 
dow e{E) vanishes includes a growing number 0(ff 1 / 4 ) of modes. The result 
(with equivalent choices of L(E) and e(E)) has been proved for Schrodinger 
operators with smooth confining potential using coherent states [ 22 ] ■ Con¬ 
jecture 0] is a stronger version of the spectral measure theorem m eh ] 


lim 

E —^oo 


l 

N(E) 


^2 \{(pn, Mm)\ 2 

m,n:E n £[0,E] 

Oi<fkm k n <(3 



C A {u) 

2tt 


du, 


( 12 ) 


which holds for a < (3 independent of ergoclicity. It is known Va(E\uj) van¬ 
ishes for ergodic weak-mixing flows (weak-mixing ensures Ca (oj) is a bounded 
function), but without any proven rate |SBJ. The most accurate previous nu¬ 
merical test of Conjecture 0 is the 10% agreement found for billiards with 
A a singular boundary operator [8]. Other tests have included the quartic 
oscillator [3], bakers map BE!, and limaqon billiard m- 
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Existing numerical studies of all the above conjectures share the features 
of low and poorly-quantified accuracy (he. lack of statistical rigor in the 
tests), and relatively low mode numbers (n 10 3 to 10 4 ). In this work 
we remedy both these flaws by performing a large-scale study using non¬ 
standard cutting-edge numerical techniques which excel at very high eigen¬ 
values. In Section 12.21 we review heuristic arguments for Conjectures [H [21 
and 01 Then in Section 01 (which refer to the three Appendices) we outline 
numerical methods; we emphasize there are several innovations. However the 
reader interested in the results on the four Conjectures, and discussion, may 
skip directly to Section 01 We summarize our conclusions in Section 01 

2 Heuristic arguments for ergodicity rate 

Here we review some arguments for Conjectures 01121 an d 01 from the physics 
literature. We feel these are appropriate since quantum chaos is somewhat 
of a cross-over area between mathematics and physics. 

2.1 Random wave model 

Berry na put forward the conjecture that chaotic eigenfunctions should look 
locally like a superposition of plane waves of fixed energy E , traveling in all 
directions, with random amplitudes and phases (compare Figs. 01 and QJ). We 
will show (following [21!) diat this model satisfies Conjecture Q with certain 
7 and a. We note that Zelditch has more rigorously considered models of 
random orthonormal bases on manifolds, and shown that they possess similar 
ergodicity properties to quantum chaotic eigenfunctions [571 EH- 

Using the notation (•) to denote averaging over the ensemble of plane 
wave coefficients, the two-point correlation of a random-wave held is |K| 

(0(ri)0(r 2 )) = yoI |^ J o(fc|ri - r 2 |), ( 13 ) 

where we use the normalization vol(f2)(0 2 ) = 1 appropriate to the billiard 
area. This applies to both real (time-reversal symmetric) and complex (non 
time-reversal symmetric) waves; from now we stick to the real case. Because 
the model is a statistical one, it is meaningful to speak of the variance of a 
particular matrix clement within the ensemble, namely var(A nm ) where we 



... .J&. %• 



Figure 3: Density plot of |0 n | 2 for an eigenfunction with /c n = 999.90598 • • •, 
that is, E n ~ 10 6 , and level number n ~ 5 x 10 4 . There are about 225 
wavelengths across the diagonal. 
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Figure 4: Density plot of one sample from the ensemble of random plane 
waves with the same wavenumber magnitude |k| = k and mean intensity as 
the eigenfunction in Fig. |2I shown in a square region of space. (Note there 
are no boundary conditions imposed). The ‘stringy’ structures visible to the 
eye are a feature of the constant wavevector magnitude; they disappear if a 
finite range of k is included nr 
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use A nm := (cf) n , A<j) m ). A simple calculation using (THU) and Wick’s theorem 
for Gaussian random variables gives 


( 0 n (ri) 0 m (ri) 0 n (r 2 ) 0 m (r 2 )) 


vol(f 2) 5 


[4im “I - Qnm Jq ( k | F i F 2 |)J 


(14) 


where 


Qnm 


9 , n = m , 
1 , n 7 ^ m. 


(15) 


Here the symmetry factor is g = 2; we will shortly see it gives the ratio 
of diagonal to off-diagonal variance. A key assumption made was that (p n 
and (j) m are statistically independent members of the ensemble when n ^ 
m, a natural one in the context of random matrix theory (RMT) (g — 2 
for time-reversal symmetry is a standard result in the Gaussian Orthogonal 
Ensemble CD- The independence assumption, like the random wave model 
itself, remains a heuristic one, albeit one with numerical support. Clearly 
orthogonality © dictates that, when considered as functions over all of 0 , 
4> n and (j) m cannot be independent! (As was mentioned there are improved 
models which take this into account EH)- An intuitive argument can be 
made if a restriction to a small subregion of 0 is made: the eigenfunctions 
behave like random orthogonal vectors, and projections of such vectors onto 
a much smaller-dimensional subspace become approximately independent. 

Now we use m to evaluate the variance diagonal and off-diagonal matrix 
elements, writing the variance as the mean square minus the square mean, 


var(A nm ) = ( 


0n(r)0 m (r)A(r)dr 


Jn 

Qnm 


) — {A nm )" 


vol(G ) 2 

Qnm 

7T VOI(G) 



J 0 (k n \ri - r 2 |) Jo(k m \ri - r 2 \)dr 1 dr 2 


n Jn 


£- 1/2 f [ A{ ri )A(T 2 ) 



nJn l r i — r 2 | 


dr\dr 2 , 


(16) 


where we used (A nm ) = S nm A. In the final step two approximations have 
been made: i) L\k n — k m \ 1 where L is the largest spatial scale of A(r), 
meaning that the two Bessel functions always remain in phase so can be set 
equal, and ii) the asymptotic form Jo(x) ~ (2/nx) 1 ^ 2 cos(a; — 7 t/4 ) was used, 
and cos 2 replaced by its average value |, giving a semiclassical expression 
valid when kl 1 , where l is the smallest relevant spatial scale in A(r). 
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Considering the diagonal and off-diagonal cases, m implies 


var(.4„„) = g var(/l„ m ) 


( 17 ) 


in a region n ~ m close enough to the diagonal. The diagonal case of CED 
gives the power law 7 = 1 / 2 , 


V a (E) S3 a „£r 1/2 , 


(18) 


where the prefactor takes the form of a Coulomb interaction energy of the 
‘charge density’ A(r), 



(19) 


Note that this model takes no account of the billiard shape or boundary 
conditions. 

2.2 Classical autocorrelation argument (FP) 

Fcingold and Peres ra were the first to derive a semiclassical expression for 
diagonal variance in chaotic systems. There are two steps: 

(i) relating off-diagonal variance to Ca(oj) the power spectral density, yield¬ 
ing Conjecture IU then 

(ii) relating diagonal variance to off-diagonal variance close to the diagonal, 
yielding Conjecture El 

Step (i) : Our presentation is loosely based on Cohen [2H] : we emphasize 
that this is not a mathematical proof, rather one form of a heuristic common 
in physics literature (cf. jSH ES GU 1211). The autocorrelation m of the 
‘signal’ A(r(t)) associated with a uniformly-distributed ergodic unit-speed 
trajectory r(i) is Ca(t) = A( 0 )A(r), where the ergodic theorem was used 
to rewrite the time average as an average over initial phase space locations 
(r o ,0 o )- Fixing r gives the function of phase space -4(r o ,0 o ) := A(0)A(r). 
Applying QET to this function A, gives as E —> 00 , 





( 20 ) 


n:E n £[E,E+L(E)] 
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where A(t) is the quantization of A, shifted in time (according to the Heisen¬ 
berg picture of quantum mechanics). A window L(E ) = 0(E 1 ^ 2 ) is sufficient 
for validity of QET |63| . Assuming a wave dispersion relation u = k the 
operator A(t) is expressed in the eigenfunction basis, 

(0 n ,i(t)0 m ) = A nm e- i{km ~ kn)t . (21) 

Using this and inserting a sum over projections onto all eigenfunctions into 
( 1201 ) gives 


Nl(E)- 1 \Anm\ 2 e~ i{ - krn ~ kn)r ~ C a (t). ( 22 ) 

n:E„e[E,E+L(E)] m= 1 

Taking the inverse Fourier transform of the definition in m, recognising 
lim e ^ 0 (2e) _1 l[-e,e] as the Dirac delta function, and using (1221) . gives 

^ J e~ iuT V A {E] u) du ~ ^lc A (r). (23) 

Finally, taking a Fourier transform and substituting A k{E) gives Conjec¬ 
ture 01 

Step (ii): To approach the diagonal we take the limit u —> 0; here 
C a {uj) is well-dehned and bounded since the flow is weak-mixing [60|. The 
expectation that diagonal variance should exceed off-diagonal variance by 
the time-reversal invariance symmetry factor (El) with g — 2 gives Conjec¬ 
ture [2J Feingold -Peres [213 justify this by considering <f)± := (0 n ± </> m )/\/2 
for eigenfunctions cf) n and (j) m with sufficiently small E n — E m . Matrix el¬ 
ements are then (c f> +1 A(j)_) = ( A nn — A mrn )/ 2, where A = A* was used; 
these are expected (by rotational invariance) to have the same variance as 
A mn . Treating these quantities as independent statistical variables, taking 
variances, and recognising that mean values A cancel on the right, gives 
var(A mn ) = var(A + _) = var(A nn )/2, that is, g — 2. Although this was 
not discussed by FP, the system must be assumed to be without further 
symmetry, or other values of g may result. 

Finally, the physical significance of the correlation functions used above 
should be noted: the right-hand side of m is proportional to dissipation 
(heating) rate in a classical system driven at frequency u with the forcing 
function A, and the left-hand side to quantum dissipation rate under equiva¬ 
lent forcing (with A) within linear response theory (for reviews see [20, Hi!?]). 
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Figure 5: Demonstration that none of the 6812 inodes are missing in the 
wavenumber window k n G [650, 750]. The level counting function © is plot¬ 
ted after the first two Weyl terms N Weyl (E ) = (vol(Q) / An) E — (\Y\/A-k)^/E 
have been subtracted. The horizontal axis shows wavenumber k = E 1//2 . 
Spectral rigidity ensures that a single missing or extra mode can be de¬ 
tected |2B]; this would be visible as a permanent jump of size 1 (the gap 
between the dotted horizontal lines). No such jump occurs. 
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3 Computing quantum matrix elements 

Our nonsymmetric billiard 0 is shown in Fig. [TJi and defined in the caption. 
Its classical autocorrelation has already been studied jjj]. In App. [3 we 
present the method used to compute Ca{oj) to an accuracy of a fraction 
of a percent (this is much less time-consuming than the following quantum 
computations). Note that since two walls are the coordinate axes it is a 
desymmetrized version of a billiard of four times the area formed by unfolding 
reflections in x and y. Desired eigenfunctions are then the subset of odd-odd 
symmetric eigenfunctions of the unfolded billiard, a fact that enables use of a 
smaller symmetrized basis set (hence higher eigenvalues), and a reduction of 
the effective boundary to r desym , the part of Y which excludes the symmetry 
lines (axes), accelerating the quantum calculation by a factor of about 4 [7j. 

We chose the test function A(r) to be the characteristic function of the 
subdomain 14^ shown in Fig. which falls one side of the straight line 
<914^\T. Our choice of the shape of 14^ was informed by the issue of boundary 
effects raised by Backer et al. [5j, the main point being that within a boundary 
layer of order a wavelength, there are Gibbs-type phenomena associated with 
spectral projections, and by choosing a large angle of intersection of the 
line with Y their contribution is minimized. Our classical mean is A = 
vo1(14a)/ vol(14) ~ 0.55000. Matrix elements (</> n , A(p m ) = (</> n , (p m )n A are 
computed using integrals of eigenfunctions in an efficent manner described 
shortly. 

Eigenfunctions and eigenvalues were found with the ‘scaling method’ Ik>I 
IT., outlined further in App. El This is a little-known basis approximation 
method, a variant of the Method of Particular Solutions 0021 , which uses 
collocation on r desym to extract all eigenvalues lying in an energy window 
E n G [E,E + L(E)\, where L(E) is 0(E ly/2 ). The equivalent wavenumber 
window is 0(1). The spectrum in larger intervals can then be found by 
collecting from sufficiently many windows. We are certain that all modes 
and no duplicates have been found in the desired intervals, as demonstrated 
by comparing the counting function against Weyl’s law in Fig. E At energy 
E (wavenumber k ) the required basis size N = 0(N SC ), the ‘semiclassical 
basis size’, is 

N sc := (24) 

7T 

where |r desym | indicates the desymmetrized perimeter. Rather than evaluat¬ 
ing (j) n at a set of points covering all of 14, which would be very expensive, only 
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basis coefficients of (j) n are computed, from which f> n (r) can be later computed 
at any desired location. Computational effort is 0(N 2 ) = 0(E) per mode, 
assuming 0 (N 3 ) theoretical effort for dense matrix diagonalization (in fact 
because of memeory limitations, for large N it was not this favorable). Be¬ 
cause of the simultaneous computation of many modes, the scaling method 
is faster by 0 (N ) than any other known method (see overview in {221 7]) 
such as boundary integral equations B]. We needed N ~ 3500 at the largest 
wavenumber reached (k n cs 4000, at n ~ 7 x 10 5 ) for the billiard under study; 
in this case the resulting efficiency gain is roughly a factor of a thousand! 
Similar efficiency gains have been reported in other studies of the Dirichlet 
eigenproblem at extremely high energy in both 2D 23 1E2> 13 and 3D j^U] . 
Only a couple of studies in billiards have computed eigenfunctions at greater 
n, and they invariably involved shapes without corners (for example TT )| ). 
Note that in App. El we outline the basis set innovation that allows us to 
handle non-convex shapes and (non-reentrant) corners effectively. Despite 
its success the scaling method has not yet been analysed in a rigorous fash¬ 
ion JITT . 

Once a large set of eigenfunctions (in the form of their basis coefficients) 
have been found, such as the example plotted in Fig. E matrix elements may 
be efficiently computed as we now show. The following is proved in App. IQ 

Lemma 3.1 Fix E > 0 and let —A u = Eu and —Av = Ev hold in a 
Lipschitz domain VLa £ M 2 - Let n be the outwards unit normal vector at 
boundary location r e 8 Qa ( r is measured relative to some fixed origin), and 
ds be surface measure. Then 



{u,v)q a = —— ® (r-n)(Euv — Vu-Vn) + (r-V«)(n-Vn) + (r-Vn)(n-Vit) ds. 


This boundary integral identity, with the substitution u = v = <f n , allows 
diagonal matrix elements to be calculated using ID rather than 2D numerical 
integration. Since typically 10 quadrature points per wavelength are needed 
for integration, and our system is up to hundreds of wavelengths in size, this 
is an enormous efficiency gain of order 0(N) or 10 3 . Note that the boundary 
integrand is nonzero only on the line 8 Qa \ T. Off-diagonal matrix elements 
are found via the following identity which is a simple consequence of the 
Divergence Theorem. 
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Best-fit 

a (1/2) 

u bf 

Random wave 

®RW 

Feingold-Peres 

n F p 

prefactor 

0.334 ±0.003 

0.5995 ±0.001 

0.3550 ± 0.0004 

deviation from aB F /2) 

— 

79 ± 1% 

6.5 ± 0.9% 


Table 1 : Comparison of best-fit prefactor a in Conjecture [T] assuming fixed 
power 7 = 1/2, against the random wave prediction (THfll . and the Feingold- 
Peres prefactor a FP := 2CA(0)/ vol(fl) of Conjecture [21 

Lemma 3.2 Let — A u = E u u and —Av = E v v hold with E u ^ E v , and other 
conditions as above, then 

( u , v)n A = — -— ® (nn ■ Vu - vn ■ Vu) ds. 

Eu E v Jgci A 

Again we use u = (j) n , E u = E n , v = (j) m , E v = E m , and the integrand is 
nonzero only on \ T- 

Thus values and first derivatives of eigenfunctions on boundaries alone 
are sufficient to evaluate all matrix elements. Note the eigenfunctions need 
never be evaluated in the interior of 12^. For the boundary integrals, 0(N ) 
quadrature points are needed, and at each point 0(N ) basis evaluations are 
needed to find (j) n and its gradient, giving 0(N 2 ) effort per eigenfunction, the 
same effort required to fold modes by the scaling method. The calculations 
reported in this work took only a few CPU-days (1GHz Pentium 111 equiva¬ 
lent, 1-2 GB RAM) in total. The effort is divided roughly equally between 
evaluating basis (Bessel) functions and their gradients at the quadrature 
points, and dense matrix diagonalization. 

4 Results and discussion 

4.1 Diagonal variance (Conjectures [T| and [^) 

Fig. [HI shows a sample of raw diagonal matrix element data, using only eigen¬ 
values in certain intervals, up to k & 10 3 (E & 10 6 ). From this we chose 
a sequence of E values and computed Va(E ) at each. Computing matrix 
elements at high eigenvalues is very costly, so it would be inefficient and im¬ 
practical numerically to grow the interval width as L(E) = cE 1 / 2 for some 
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Figure 6: Scatter plot of diagonal matrix elements {(f> n , A(f) n ) plotted against 
energy eigenvalue E n . The gaps are due to the fact that only certain windows 
on the E axis have been computed; within each window all eigenvalues are 
found. The windows shown in the main plot correspond to wavenumbers 
k n G [100,500], k n G [650,750] and k n G [950,1050], giving a total of 28171 
modes. The classical mean A is shown as a horizontal line. The inset shows 
the complete energy range (including the 2718 higher modes not shown in the 
main plot) on a log scale, with larger points to make extreme values evident. 
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Figure 7: Log-log plot of diagonal variance Va(E) (with ±lcr errorbars) as 
a function of energy E. The smallest errorbar is about 1.4% (at E ss 10 6 ). 
For Conjecture [D the best-fit power law (7bf, a bf) is shown as the dashed 
line. The random wave prediction m and m is the dash-dotted line. 
Conjecture [21 (FP) is shown as a solid line, and the ratio of Va(E) to this 
prediction is shown in the inset. 
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constant c. This would either involve ignoring, or breaking in to very short in¬ 
tervals (which would introduce large relative fluctuations) low eigenvalue data 
which is cheap to collect, or requiring vast numbers of inodes at high eigen¬ 
value which is too expensive. Rather, we chose a convenient value of L(E) 
at each E that allowed a large number of modes M to be summed at that E, 
while still allowing access to the highest E values possible ( E & 1.6 x 10'). 
Specifically, we split the lowest interval shown (k n G [100,500]) into intervals 
containing about 10 3 modes each, kept the two other intervals shown intact 
(each containing an extremely large number of order 10 4 ), and at the higher 
eigenvalues (see inset of Fig. m chose intervals containing 200-700 modes 
each. 

One may ask how the Va{E) values obtained this way would differ from 
those obtained using a strictly growing L(E) = cE 1 ^ 2 . To indicate the flu- 
tuations in Va(E) expected from summing over a finite sample size M, we 
included an errorbar of relative size a/2 /M. For example M = 20000 corre¬ 
sponds to 1 % errorbar (illustrating that high statistical accuracy is computa¬ 
tionally intensive). This model does assume that deviations A(f> n ) — A\ 
are roughly statistically independent, an assumption motivated by RMT B2I- 
We find (see inset of Fig. □) that observed fluctuations in Va(E) fit this as¬ 
sumption well, so that even if correlations are present they do not affect 
our conclusions much. Thus if longer sums were indeed computed using 
L(E ) = cE 1 / 2 , it is likely that their values would lie within the errorbars 
shown. The resulting local variance is presented in Fig. 0 A power-law 
dependence is immediately clear. We fitted the power-law model of Conjec¬ 
ture [TJ obtaining a best-fit power 

7 BF = 0.479 ±0.009, (25) 

differing by only 4 ± 2 % from the random-wave model and Conjecture [21 
value 7 = 1/2. For this fit we used weighted least-squares, weighted using the 
numbers of modes M (equivalently, errorbars) for each interval. In a standard 
maximum-likelihood framework | 0 H j we marginalized over a to obtain the 
quoted errorbar in 7 . We also excluded a low-eigenvalue regime found to 
be non-asymptotic, taking only data with E > E min = 1.6 x 10 5 . (Our 
criterion here was that if any lower E min was used, the fitted 7 BF was found 
to dependent on E min )- This best-fit power-law is shown in Fig. [3 the data 
are completely consistent with Conjecture Cl 

The random wave model and Conjecture [21 both involve the power 7 = 
1/2. Therefore to test their validity the power was held fixed at this value 
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while fitting only for the prefactor a, with results given in Table |T] Notice 
that the random wave model is a poor predictor; this might be expected since 
this model takes no account of the boundary conditions, yet the support of 
the test function A extends to the boundary and covers a large fraction of 
the volume. Intuitively speaking, images (boundary reflections) of A are not 
being taken into account, and evidently this is a large effect. 

The prefactor (which is dominated by the large numbers of modes at 
E « 10 6 ) is overestimated by Conjecture [21 by only 6.5%, however this is 
statistically significant (a 7cr effect). The systematic deviation is highlighted 
in the inset of Fig. [71 where it appears that the deviation decreases with 
E. This is especially clear for E ps 10 6 or below, and asymptotic agreement 
with Conjecture [21 is not inconsistent with the large scatter of the highest six 
datapoints (since they have larger errorbars). If Conjecture [21 is correct, our 
results show that convergence must be alarmingly slow. This may explain 
why previous numerical billiard studies II El found various power-laws and 
prefactors differing by up to 20% from Conjecture [21 depending on choice of 
billiard and A: they had failed to reach the asymptotic regime. We achieve 
this with mode numbers 100 times higher than these studies. We might 
model slow convergence to Conjecture [21 by a correction, for example, 

V A {E) ~ l 1 - bE ~ 0 + °(^)) > (26) 

with (3 sufficiently small, and b sufficiently large. Our current data do not 
allow a meaningful fit for (3 and 6, however, they suggest 0.25 < f3 < 0.5, and 
that b is several times greater than 1. We note that in the physics literature 
periodic-orbit corrections to Conjecture |21 derived by Wilkinson M were 
later claimed by Prosen not to contribute at this order |39|. Clearly more 
numerical and theoretical work is needed on this issue of slow convergence. 

4.2 Quantum Unique Ergodicity (Conjecture [31) 

We now examine individual diagonal matrix elements. The inset of Fig. |o] 
enables extreme values to be seen: it is clear that there are no anomalous 
extreme values which fall outside of a distribution which is condensing to the 
classical mean. Both the maximum 0.6811 and minimum 0.3437 of ((p n , A<fi n ) 
occur at E n < 2 x 10 4 , visible at the far left side. This is strong evidence for 
QUE in this system. Since about 30000 modes are tested, the density of any 
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Figure 8: Histogram of rescaled diagonal deviations z n := A(ft n ) — 

H)/Va(-En) 1//2 , using the best-fit power law form for Conjecture [TJ compared 
against an appropriately-normalized Gaussian distribution of unit variance. 
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excluded sequence can be given an approximate upper bound of 3 x 10 5 . We 
note that it is possible (althouth unlikely) that by an unfortunate choice of A, 
non-uniform or scarred modes occur which do not have anomalous (</> n , A(j) n ) 
values. The only way to eliminate this possibility would be to repeat the 
experiment with a selection of different A functions. 

What is the distribution that the deviations (cp n , A(p n ) — A follow? We 
have rescaled these deviations according to the best-fit form of the variance 
in Conjecture |TJ and histogram the results in Fig. El The distribution is 
consisitent with a gaussian, with an excellent quality of fit. (The slightly 
fatter tail on the low side is entirely due to low-lying modes; as can be seen 
in Fig. El these have a skewed distribution). 

4.3 Off-diagonal variance (Conjecture HI) 

We decided, for reasons of numerical practicality, to test Conjecture El over 
a range of to, but at a single (large) value of E. This was performed using 
the single sequence of 6812 eigenfunctions with k n E [650,750], from which 
about 2.3 million individual off-diagonal elements (</> n , A(j) m ) were calculated, 
namely those in the block n, m : k n , k m E [650, 660], the block n, m : k n , k m E 
[660, 670], etc, with 10 blocks in total. The first such matrix block is plotted 
in Fig. El note the strong diagonal band structure (‘band profile’ |20|). The 
mean off-diagonal element variances lying in successive intervals of width 
0.1 were then collected. Thus we are testing Conjecture El with a window 
of L{E) ~ 1.4 x 10 5 and e(E) = 0.05. Onr extremely large choice of L(E) 
allowed statistical fluctuations to be minimized (errorbars were estimated as 
in Section EH and are generally less than 1%). We believe this is the most 
accurate test of the conjecture ever performed. 

The resulting band profile Va{E\uj) is compared against Conjecture El hr 
Fig, mu The agreement is excellent, with generally less than 3% discrepancy. 
Let us emphasise that there are no fitted parameters. There appears to be 
statistically significant deviations: the peaks and valleys (points of highest 
curvature, including uj — 0) are exaggerated more in the quantum variance 
than the classical. Care has been taken to ensure that this was not due to 
numerical errors in Ca(oj) (it was smoothed only on a much finer scale, see 
Appendix \M As with the diagonal variance, this may reflect slow conver¬ 
gence. The diagonal quantum variance, divided by g = 2, is also shown. If 
step (ii) of Section 12721 applied exactly then this would coincide with Va(E , 0), 
however a 3± 1% difference is found. It is unlikely but possible that this 
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Figure 9: Density plot of block of squared matrix elements \{(p n , A(p m ) — 
5 nm A\ 2 for the 637 modes lying in k G [650,660]. The range white to black 
indicates zero to 1.7 x 10~ 3 . Individual elements appear uncorrelated, the 
only visible structure being the intensity varying with spectral measure (band 
profile). 
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Figure 10: a) Off-diagonal quantum variance Va{E',u> ) (errorbars) at a single 
E ss 4.9 x 10 5 , as a function of distance oj (wavenumber units) from the 
diagonal, compared against Conjecture 0 estimated using the classical power 
spectrum (solid line). Errorbars are 0.7% for quantum (near the diagonal), 
0.2% for classical. Also half the diagonal variance is shown as two horizontal 
dotted lines indicating ±1 standard deviation. Inset is a zoom on the u —> 0 
region, b) Relative deviation (difference from 1 of ratio) between quantum 
variance and Conjecture |3 
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is merely a statistical fluctuation (a null result). Thus, at this energy, the 
diagonal variance prefactor discrepancy of 7% seems to result from the ad¬ 
dition of two roughly equal effects: step (i) has about 4% discrepancy, and 
step (ii) about 3%. Our data does not rule out the possibility that g is closer 
to 1.85 than to 2, which would imply slight positive correlations between 
neighboring eigenfunctions. 

In terms of individual elements, we find no anomalously large values. 
This strongly supports ‘off-diagonal QUE’: the vanishing of every single off- 
diagonal element, a stronger result than i;52|. A preliminary examination 
suggests uncorrelated gaussian distribution of elements, with variance given 
by cn-distance from the diagonal, but we postpone analysis for future work. 

4.4 Discussion on existence of scars 

In the physics community the existence of scars is well-known, as are the¬ 
oretical (non-rigorous) models. Heller p])) EDj put forward a semiclassical 
explanation based on enhanced short-time return probability for wavepack- 
ets launched along the least unstable periodic orbits (UPOs), which has been 
elaborated na na in E3 - Although the meaning of ‘scar’ varied historically, 
it is now taken to mean any deviation from the random wave prediction of 
eigenfunction intensity near a UPO EH- Scar ‘strength’ depends on what 
test function you use to measure it R Bj : in physics this test function is com¬ 
monly not held fixed as the limit E —> oo is taken, rather it is chosen to 
collapse microlocally onto a UPO with a coordinate-space width ~ EU 1 / 4 . 
By this measure, scar strength is believed not to die out in the semiclassical 
limit Typical scar intensities (|</> n | J along the UPO) do not decay, but 
their width, hence the associated probability mass, vanishes as 0(EU 4 / 4 ). 

However, in the mathematics community questions of uniformity of eigen¬ 
functions and QUE are presented in the form of weak limits, that is, limits of 
matrix elements of fixed, O th -order pseudo-differential operators. Persistence 
of scarring is taken to mean existence of a subsequence with | ( (p n , A(p n ) —A\ = 
0(1). Thus we might distinguish a physicist’s scar , where a probability mass 
vanishing as 0(EU 1//4 ) is associated with the UPO, from a mathematician’s 
scar (or ‘strong scar’) which carries 0(1) probability mass (as in SI]). 

Our results (Section 14.11) suggest strong scars do not persist asymptoti¬ 
cally, but is consistent with the persistence of physicist’s scars, in fact giving 
the same power-law |(0 n , Acj) n ) — A | = 0(EU 1//4 ) expected from scar width. 
Heller’s numerical demonstrations of apparently strong scarring were done 
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at mode numbers n ~ 2 x 10 3 , which, in light of our work, is well below 
the asymptotic regime. It is now believed by physicists that in 2D Anosov 
billiards strong scarring does not persist [U (S3 ESI, but there still exist con¬ 
troversies about the width of scars [ 23 ], and in related quantum models the 
mechanism of scarring is an active research area [^J. Mathematically, the 
issue remains open. 


5 Conclusions 

We have studied a generic Euclidean Dirichlet billiard whose classical dy¬ 
namics is Anosov, and studied quantum ergodicity of both diagonal and 
off-diagonal matrix elements. By accessing very high eigenvalues (100 times 
higher than previous studies) using the scaling method and boundary inte¬ 
gral formulae for A(p m ) in the case of piecewise-constant test function 
A, we believe we have reached the asymptotic regime for first time. We also 
have unprecedented statistical accuracy due to the large number of modes 
computed. A summary of evidence found for the four conjectures from the 
Introduction is as follows: 

• Conjecture |U : Diagonal variance shows excellent agreement, with 
power-law 7 = 0.479. 

• Conjecture |21 : Diagonal variance is consistent with the Feingold-Peres 
prediction, but with quite slow convergence (for instance, a 7% over¬ 
estimate remains at E — 10 6 ). This contrasts a random-wave model, 
which overestimates the fitted prefactor by 80%. 

• Conjecture 01 : Compelling evidence for QUE in this system (density of 
exceptional subsequence < 3 x 1CT 5 ). 

• Conjecture 0] : Excellent agreement for off-diagonal variance (of order 
3% discrepancies at E = 5 x 10 5 ), and evidence for off-diagonal QUE. 

Performing this large-scale study required limiting ourselves to one bil¬ 
liard D and one operator A. In order to complete the picture, the discrepancy 
between the predictions of Conjectures [21 and 0 and numerically-measured 
variances should be studied as a function of f2 and A. 
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A Classical power spectrum 

We use standard techniques 1271 to estimate CAM- For a particular trajec¬ 
tory, launched with certain initial location in phase space, A(t) is a noisy 
function (stochastic stationary process). We define its windowed Fourier 
transform as 

A(u) := [ A{t)e iujt dt , (27) 

Jo 

where the window is a ‘top-hat’ function from 0 to T. Using with (HDD and ©, 
and taking care with order of limits, we have the Wiener-Khinchin Theorem, 

CUM = Jim A*(uj)A(uj). (28) 

T —too Z7T1 

For this single trajectory, A(u) is a rapidly-fluctuating random function of 
u, with zero mean (for uj A 0), variance given by 2 / kTCa {M and correla¬ 
tion length in uj of order 27r/T. (As T —* oo, the ^-correlation becomes a 
delta-function). Thus (j2HD converges only in the weak sense, that is, when 
smoothed in uj by a finite width test function. 

A given trajectory is found by solving the particle collisions with the 
straight and circular sections of T, and A(t) is sampled at intervals At = 0.02 
along the trajectory (recall we assume the particle has unit speed). Then 
A(u) is estimated using the Discrete Fourier Transform (implemented by an 
FFT library) of this sequence of samples, giving samples of the spectrum at 
uj values separated by Aoj = 2n/T. The correlation in uj is such that each 
sample is (nearly) independent. At was chosen sufficiently short that aliasing 
(reflection of high-frequency components into apparently low frequencies) 
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was insignificant. A trajectory length T = 10 4 (about 1.8 x 10 4 collisions) 
was used. The finiteness of T causes relative errors of order f corr /T, where 
t COIr ~ 2 (for our domain) is the timescale for exponential (since the billiard 
is Anosov) decay of correlations. Thus more sophisticated window functions 
are not needed. 

Given A{u) we use (0H1) . with the T given above, to estimate Ca{oj). We 
smooth in a; by a Gaussian of width u; sm = 0.03. This width is chosen to be as 
large as possible to average the largest number of independent samples from 
the neighborhood of each u, but small enough to cause negligible convolution 
of the sharpest features of C-a{oj). 

Finally, in order to reduce further the random fluctuations in the estimate, 
n r = 6000 independent trajectory realizations with random initial phase 
space locations were averaged. An estimate for the resultant relative error 
e in Ca{w) can be made by counting the number N of independent random 
samples which get averaged, and using the fact that the variance of the square 
of a Gaussian zero-mean random variable (i.e. x 2 distribution with 1 degree 
of freedom) is twice the mean. This gives 
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which numerically has been found to be a conservative estimate. In our case 
e ~ 2 x 10” 3 , that is, about 0.2% error. 

The zero-frequency limit Ca( 0) is found using the smoothed Ca(uj) graph 
at to = 0, and therefore is an average of frequencies within 0(u; sm ) of zero. 
This is justified because correlation decay (weak mixing) causes all moments 
of Ca(t) to be finite, hence there is no singularity in Ca(w) at u = 0 (it can 
be expanded in an even Taylor series about u = 0 with finite coefficients). 


B Scaling method for the Dirichlet eigenprob- 
lem 

The scaling method for the solution of the Dirichlet eigenproblcm in star¬ 
shaped domains was invented by Vergini and Saraceno ESj, and considering 
its great efficiency it has received remarkably little attention. Here we give 
only an outline. 
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The method relies on the remarkable fact that the normal derivatives of 
eigenfunctions lying close in energy are ‘quasi-orthogonal’ (nearly orthogonal) 
on the boundary, with respect to the boundary weight function r n := r ■ n. 

Lemma B.l Let G W l , d > 2 be a Lipshitz domain with boundary T and 
Dirichlet spectral data {£)}, {<f>j}. Let n and r be defined as in Lemma 1X71 
with r n := r • n. Then, for all i > 1, j > 1, 

Qij ■= j> r n (n ■ V 0 i)(n ■ ds = 2E i 8 ij + ^ E% ^ r 2 0 j ) Q) ( 30 ) 

This is proved in App.P A corollary is that, since r 2 := r-r is a bounded op¬ 
erator on the domain, off-diagonal elements of Q must vanish quadratically 
as one approaches the diagonal. Thus the matrix with elements Qij/2Ei 
approximates the identity matrix, when restricted to an energy window 
E tl Ej G [E — Co, E + eo], if the window size remains relatively narrow 
e 0 = o(E l/2 ). 

We choose a ‘center’ wavenumber k = E 1 / 2 , near which we are interested 
in extracting eigenfunctions, and relative to which the wavenumber shift of 
mode i is ufik) := k — ki. Consider an eigenfunction </>* for which uy < 0 
and \uf -C 0(1). We create a version spatially rescaled (dilated about the 
origin) by an amount needed to bring its wavenumber to k, that is, y,:(fc, r) : = 
(pfikr/ki). We call this function fc-rescaled. Thus we have —A Xi = E\i 
everywhere inside fl, with r) = 0 on the rescaled boundary (he. for all 
kr/ki G T). The rescaled eigenfunction can be Taylor expanded in u 

Xi(k, r) = + = ^(r) + ■ V0i + 0(o; 2 ) 

= jpr n n ■ + 0(u> 2 ) for r G T, (31) 

where Dirichlet boundary conditions were applied. We construct a basis of 
N functions ffik, r), satisfying —A£; = Efj inside no particular boundary 
conditions on T, and non-orthogonal over fl We assume they approximately 
span the linear space in which rescaled eigenfunctions live, so that 

N 

Xi(k, r) = E Xufifik, r) + 6j(r) for all i of interest, (32) 

1=1 

where the error e t can be made negligibly small for some N. In practise 
N need exceed N sc defined in (El by only a small factor (2 or less). Our 
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goal is then to solve for a shift ux and the corresponding i th column of the 
coefficient matrix X. We can do this by simultaneous diagonalization of 
quadratic forms. We define two symmetric bilinear forms on the boundary, 

® — uvds, (33) 

Jr r n 

\ ® — (ur ■'Vv + vr ■ Va) ds. (34) 

k J r r n 

Note that defining these forms brings the extra requirement that the domain 
be strictly star-shaped about the origin (r n > 0), which from now on we 
assume. In the rescaled eigenbasis / is, via m and flUl) 

f{XhXj ) — T~j~Qii + ^( tc ’ 3 ) = + 0(u; 3 ), (35) 

K'i rZj 

a matrix which is close to diagonal, because of the closeness of Q to the iden¬ 
tity. In the same basis, recognizing that for fc-rescaled functions g is equiva¬ 
lent to df /dk, the derivative of (1341) with respect to the center wavenumber, 
and using duJi/dk = 1, we have 

g(Xi, Xj) = — 7 * J Qij + 0{^~) = 4:UJi5ij + 0(u> 2 ), (36) 

rCi rCj 

so g is also close to diagonal. Thus the set {xj} with small uy | approximately 
diagonalizes both bilinear forms, with the approximation error growing as a 
power of \u>i\. As we explain below, in practise the converse applies, that is, 
by simultaneously diagonalizing / and g we can extract the set of eigenfunc¬ 
tions lx*} with smallest \uii\. Therefore, loosely speaking, when the boundary 
weight function 1 jr n is used, domain eigenfunctions are given by the simul¬ 
taneous eigenfunctions of the (squared) boundary norm and its /c-derivative. 

We perform the diagonalization in the basis (1321) . That is, matrices 
Firn ■= /(6,^rn) and Gim \= g(Ci,Cm), with l,m = 1 ■ ■ ■ N, are filled. This 
requires basis and first derivative evaluations on the boundary. It is an el¬ 
ementary fact that given a positive matrix F and a symmetric matrix G 
there exists a square matrix Y and a diagonal matrix D := diag{/q} which 
satisfy Y T FY = / and Y T GY = I). The matrices Y and D can be found by 
standard numerical diagonalization algorithms in 0(N 3 ) time. If (1351) and 
© held without error terms, and the number of modes i for which they 
held were equal to (or exceeded) the basis size N, then we would be able 


f(u,v) := 
g(u,v) := 


31 




directly to equate the columns of Y with the desired columns of A" (barring 
permutations). In this case uy = 2//q would also hold, from which the de¬ 
sired wavenumbers ki follow. However, using Weyl’s law and (OH) is follows 
that such a large number of modes requires that the largest \uj t \ is of or¬ 
der unity, in which case errors in (EH) and EHD would become unacceptable. 
It is an empirical observation found through numerical study that in fact 
columns of Y corresponding to the largest magnitude generalized eigenvalues 
Hi (and therefore the smallest shifts |oy|) do accurately match columns of X. 
Thus perturbations by other vectors in the span of basis functions are small. 
Further discussion is postponed to a future publication d 

We mention a couple of other implementation details. Because the gener¬ 
alized eigenproblcm turns out to be singular it is truncated to its non-singular 
part mm If columns of Y are normalized such that Y 7 FY = / holds then 
the resulting eigenfunctions can be normalized over by dividing the i th 
column of Y by \f2ui. Depending on the choice of basis, spurious solutions 
can result; they are easily identified because their norm over fi, computed 
by the following Rellich identity, is not close to 1. 

Lemma B.2 (Rellich) With the definitions of Lem,m,a \B . 11 for all j > 1, 

^yV n (n ■ Vfj) 2 ds = 1 (37) 

This identity is a special case of Lemma ETT1 found by substituting Q 4 = Q. 
u — v — (j)j, and E = Ej, and recognising r • Vfj = r n n- Vfj. The maximum 
|uy| in which levels of useful accuracy are found is of order 0.2 /R where R is 
the largest radius of the domain. The lack of missing modes obtained with 
this method is illustrated by Fig. 0 There are several implementation issues 
and improvements that we do not have space to discuss here mmm- 
A word about the basis set choice {£(£;, r)}, 1 < i < N, is needed. Until 
now plane waves (including evanescent plane waves [221) or regular Bessel 
functions da have been used. These fail for non-convex domain shapes, or 
those with corners, thus to tackle the domain in this study a basis of irregular 
Bessel (he. Neumann) functions, placed at equal intervals along a curve T + 
exterior to D, was developed by the author. T + is defined as the set of points 
whose nearest distance to T is T, with kD = 7 (roughly one wavelength 
distant). This was found to handle (non-reentrant) corners successfully. It 
performs extremely well for all shapes that have been attempted so far. The 
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basis size N is about 1.5iV sc (see (El), thus, depending on required accuracy, 
about N/ 20 useful modes are found per dense matrix diagonalization (0(iV 3 ) 
effort). This is O(N) faster than other boundary methods; we remind the 
reader that N is larger than 10 3 in our work. 

C Identities involving eigenfunctions of the 
Laplacian 

Let there be constants E u > 0, E v > 0, and let —A u = E u u and —Av = E v v 
hold in a Lipschitz domain G K d , for some general dimension d > 2. The 
following expressions for the divergence of certain vector fields result from 
elementary calculus. By VVu we mean the second derivative tensor (dyad). 


V ■ (vVu) = —E u uv + Vu ■ Vu (38) 

V ■ (ruu) = duv + ur ■ Vu + ur ■ Vu (39) 

V • (rV« • Vu) = dVu ■ Vu + Vu ■ VVu • r + Vu ■ VVu • r (40) 

V ■ [(r • Vu)Vu] = Vu ■ Vu — E v v r • Vu + Vu • VVu • r (41) 

V • (r 2 uVu) = 2vr ■ Vu — E u r 2 uv + r 2 Vu • Vu (42) 


First we prove Lemma 13.11 (also see App. H of |T:). Consider the following 
four equations: flUf (SDH, (EH, and its counterpart obtained by swapping u 
and v. Integrating each of these over Qa, then applying the Divergence 
Theorem, gives four expressions for surface integrals in terms of domain 
integrals. Substituting these into the right-hand side of the expression in 
Lemma urn and setting E u = E v = E, gives, after cancellation, the left- 
hand side. 

A similar but more complicated technique proves Le mm a lB.ll First we 
consider modes i,j which are non-degenerate, that is, Ei ^ Ej. Eight equa¬ 
tions are needed: the four mentioned above, then (ETHI) and (EH an d their 
counterparts swapping u and v. Each should be integrated over and the 
Divergence Theorem applied. The following identity may then be verified by 
substitution of the eight resulting equations into its right-hand side. Using 
the abbreviations e := E u — E v , £ := E u + E v , u n = n ■ Vu, u r = r • Vu, the 
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identity to be checked is, 


.2 


/ 

Jn 



4 


u n v - v n u ) 



The substitutions u = (pi, r? = (pj , = E t and E v = Ej, and applying 

Dirichlet boundary conditions, turns the right-hand side into Qij. More 
details about how such identities are found using a symbolic matrix method 
will be postponed to a future publication na. 

Finally, the other possibility is that Ei = Ej (which need not imply 
i = j ). We take Lemma Id.II with the choices = fl, u = (pi, v = (pj, 
Ei = Ej = E, apply Dirichlet boundary conditions, and use orthonormality 
0 . Thus Lemma fB.il is proved for all choices of i,j. 
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